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ABSTRACT 

The optimized design of an exhaust emission system in 
terms of performance, cost, packaging, and engine con- 
trol strategy will be a key part of competitively meeting 
future more stringent emission standards. Extensive use 
of vehicle experiments to evaluate design system 
tradeoffs is far too time consuming and expensive. Imper- 
ative to successfully meeting the challenges of future 
emission regulations and cost constraints is the develop- 
ment of an exhaust system simulation model which offers 
the ability to sort through major design alternatives 
quickly while assisting in the interpretation of experimen- 
tal data. 

Previously, detailed catalyst models have been devel- 
oped which require the specification of intricate kinetic 
mechanisms to determine overall catalyst performance. 
While yielding extremely valuable results, these models 
use complex numerical algorithms to solve multiple par- 
tial differential equations which are time consuming and 
occasionally numerically unstable. Furthermore, keeping 
these kinetic data updated with improvements in catalyst 
formulations is also a very time consuming process. To 
resolve these issues, a simple model (SIMTWC, which 
stands for SIMple Three Way Catalyst) to predict tailpipe 
exhaust emissions has been developed by combining 
fundamental conservation expressions for mass and 
energy with a first order oxygen storage model and an 
empirical data base of steady-state catalyst performance 
obtained during engine dynamometer assessments, thus 
eliminating the need to specify the detailed kinetic mech- 
anisms. Advantages of SIMTWC over more complicated 
techniques are computational expedience and robust- 
ness with comparable accuracy. Updates to SIMTWC as 
catalyst formulations change are also accomplished in a 
timelier manner. 

This paper describes the mathematical formulation of 
SIMTWC, solution technique, input data requirements, 
and provides a comparison of model results with vehicle 
experiments. 



INTRODUCTION 

Historically, models of emission system performance 
required the specification of detailed kinetic mechanisms 
through the catalyst element to predict their performance 
[see, Oh.S. (1980), Otto,N. (1980), Montreuil, et al. 
(1992), Koltsakis, et al. (1997)]. During a calculation, 
these models specify and solve a set of partial differential 
equations to calculate the temperature and species evo- 
lution within the catalyst. The numerical algorithms are 
complex and time consuming. Herein, we take a simpler 
approach while retaining many of the key features of the 
more complex models. A first principle transient energy 
balance model is used to account for the substrate heat 
transfer dynamics, but we replace the integration of the 
species conservation equations with a steady state map 
of catalyst performance obtained during assessment on 
an engine dynamometer or reactor. Because of this fea- 
ture, the numerical integration is reduced to the solution 
of the transient energy equations and an integration to 
specify the time-dependent 0 2 storage over the catalyst 
[Gandhi, H.S. (1976)]. This results in a more robust cal- 
culation and one in which a calculation can be performed 
rapidly with accuracy similar to that produced by the 
more complex models. 

The overall objectives of this modeling approach are: 

1. To allow rapid and early design assessments of new 
concepts proposed for production exhaust/emission 
systems; 

2. To support the interpretation of experimental data 
and to serve as a means of establishing a corporate 
memory of lessons learned during experimentation; 
and 

3. To allow the easy integration of component models 
into larger system assessment models. 



Here, a "SIMPLE WORKING MODEL" approach is used 
in the development of SIMTWC. The model is simple 
enough in formulation to allow rapid calculation but 
retains sufficient detail so that the results are meaningful 
for emission system development. 
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The basic assumptions that allow this simpler approach 
are: 

1. The steady state chemistry on the catalyst is fast 
when compared to the heat transfer dynamics and 
the 0 2 storage chemistry. This quasi-steady state 
assumption implies that the reaction rates of relevant 
emission species (HC, CO, and NOx) can be repre- 
sented at any instant by an experimental steady- 
state performance map generated over a variety of 
key conditions (viz., substrate temperature, space 
velocity, air fuel ratio, and inlet gas species concen- 
trations); and 

2. Oxygen storage kinetics on the catalyst can be repre- 
sented as the filling or emptying of surface storage 
sites in correspondence with a first order kinetic 
mechanism. Since this feature happens over time 
scales that are typically longer than those associated 
with the steady state kinetics, a site balance expres- 
sion is integrated in time as well. 

Each catalyst of a given precious metal (PM) formulation/ 
loading and aging history will have a unique performance 
characteristic (or "potency") represented by its own set of 
mapping experiments. The number of experiments 
needed to characterize a particular catalyst is small 
enough to permit the rapid generation of new maps as 
catalyst formulations change/improve. 

In the SIMTWC representation, the model allows the 
determination of simple heat transfer in a pipe element 
ahead of the catalyst and the determination of tempera- 
ture and species in a catalyst element. This combination 
of elements acts as a fundamental building block to allow 
the assessment of a complete exhaust system. Gener- 
ally, these elements are combined in series/parallel com- 
binations that represent the action in a complete vehicle 
exhaust system. To use the model, precise exhaust sys- 
tem geometry, feedgas temperature, mass flow rate, A/F 
ratio, and species concentration are needed for each sec- 
ond during a typical driving cycle. 

MATHEMATICAL FORMULATION 

THE PIPE ELEMENT -The effects of the thermal mass 
and heat transfer characteristics along the exhaust sys- 
tem piping are extremely important in the determination 
of the exhaust system temperature. This is especially 
true during cold start of the vehicle when the entire sys- 
tem is warmed from ambient conditions to operational 
temperature. They are therefore important to catalyst 
light off and critical to meeting governmental emissions 
regulations, such as California LEV II or European Stage 
IV. The effects of thermal mass and convective heat 
transfer (between the gas and wall, and between the wall 
and atmosphere, plus thermal conduction along the wall) 
on the gas and wall temperature are determined using 
the conservation of energy for the gas phase (Eq. 1) and 
the surface (Eq. 2). (Note: all definitions of variables, 
units, and many physical properties are included in the 
nomenclature section). 
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where: T g , T w , and T amb are the gas, the wall and ambi- 
ent temperatures, respectively. The thermal conductivity 
of the pipe wall is A„ , and A, is its cross-sectional 

area; h t and h 0 are the internal and external heat trans- 
fer coefficients, respectively. Heat transfer between the 
exhaust gas and the inside pipe wall occurs under turbu- 
lent flow conditions and h, is determined from: 
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where: a, is an augmentation factor that can be used to 
account for increased internal heat transfer due to pulsat- 
ing exhaust flow and/or pipe bends. T f is the film tem- 
perature which equals the average of T and T w . 
D. and D n are the pipe inside and outside diame- 

' pipe r 

ters. The outside heat transfer mode in the SIMTWC cal- 
culation is approximated by an equation for free 
convection over a horizontal cylinder, and the corre- 
sponding heat transfer coefficient is: 
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a 0 is an augmentation factor to account for a small 
degree of forced convection due to air movement past the 
outside of the pipe, g is the gravitational constant, and P 

is the gas expansion coefficient and equals the reciprocal 
of absolute gas temperature. 

Boundary and initial conditions are necessary for com- 
plete problem definition and are specified as: 
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where the first boundary condition (Eq. 11) indicates that 
no heat is transferred at the ends of the system, and the 
initial condition (Eq. 13) can be used to specify an initial 
longitudinal temperature distribution along the wall. 

The above equations are used to account for the heat 
losses along the exhaust piping in the system. Currently, 
this pipe model does not include any additional chemical 
heat generation that may be generated in the exhaust 
manifold near the engine as a result of gas phase chem- 
istry under very hot conditions. However, most simulation 
results are not affected by this factor because under typi- 
cal drive cycle conditions, gas phase reaction is minimal, 
and the starting point for most simulations is downstream 
of the manifold flange beyond any significant level of 
reaction. 

THE CATALYST ELEMENT -We start our formulation 
for the catalyst with the generalized conservation equa- 
tions as seen in [Otto (1980), Li, Adamczyk and Pakko 
(1996) or Oh (1980)] and then simplify them into the form 
used in the integration of the SIMTWC model, remember- 
ing that SIMTWC is specifically designed to quickly and 
robustly analyze exhaust/ emissions systems. The follow- 
ing expressions specify the conservation of mass and 
energy in the gas phase and on the surface of the cata- 
lyst. They also represent the convective transport of spe- 
cies and energy to and from the surface as well as 
conduction along the surface, as shown in Figure 1. They 
also include chemical mechanisms to represent 0 2 stor- 
age on the catalyst surface and steady state surface 
chemistry. In the SIMTWC model, a database map of cat- 
alyst performance is used as the final representation of 
the steady state kinetics. This is done to reduce the total 
number of conservation equations to allow fast integra- 



tion of the equation set. It is also done to easily incorpo- 
rate mapping data of catalyst performance. 
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Figure 1 . Schematic of catalyst flow element 

Conservation of Energy. Species and Mass and Initial 
and Boundary Conditions: Gas & Surface 
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where e is the void volume fraction of the monolith as cal- 
culated from its geometry, v g is the gas velocity within a 

catalyst channel, and A face is the total frontal area of the 
catalyst. In the catalyst, the heat transfer between 
exhaust gas and catalyst wall is based on a laminar flow 
Reynolds number (based on the flow through a single 
catalyst channel using Eq. 20). 
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where CD is the cell density, D lm is the hydraulic diam- 
eter of a single channel in the monolith and is related to 
monolith cell density, wall thickness and cell shape 
according to the equation shown in the Nomenclature 
Section. The corresponding Nusselt number with appro- 
priate entrance effect enhancement is: 
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where z is the longitudinal distance from entrance of the 
monolith. 

To solve this equation set, one must specify several initial 
conditions for the exhaust system and one must specify 
gas inlet temperature, flow rate, and species concentra- 
tions on a second-by-second basis during the specific 
drive cycle under consideration. These conditions are 
specified as: 
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where z -zl L cat = normalized length along the catalyst. 

The Simplified Catalyst Model - We now proceed to 
reduce the equation set from this generalized form to the 
simplified form used in SIMTWC. As a first simplification 
to this equation set, the transient terms ( d I dt ) in the 
gas phase energy (Eq. 14), the gas phase species con- 
centration (Eq. 16), and the steady state surface species 
conservation (Eq. 17) expressions are neglected as in 
previous analysis [see Otto, (1980); Oh, (1980); Li, 
Adamczyk and Pakko, (1996); Koltsakis, (1997)] where 
they were shown to be small in comparison with the 
remaining terms. In addition, the definition of the void vol- 
ume fraction, e, and the geometric properties of a mono- 
lithic cell are incorporated to yield: 
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To allow rapid, robust calculation, we assume that the 
mass transfer rate is fast when compared to the surface 
reaction rate at low to mid operating temperatures; and it 
is augmented at high operating temperatures in the final 
model to accommodate conditions where the surface 
reaction is very fast and the overall catalyst conversion is 
limited by the mass transfer rate. Thus, 
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and the gas phase species (Eq. 30) and surface species 
(Eq. 31) balance equations combine to yield: 
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In order to solve Eq. 35 explicitly, one would need a 

detailed function for the reaction rate expression, R' , 
complete with rate constants, and reaction orders. How- 
ever, for simplicity and speed, we solve the species equa- 
tions in a global sense by using mapping data to specify 
the overall steady state conversion efficiency and a tran- 
sient 0 2 storage integration; 0 2 storage will be discussed 
later. The catalyst performance map generated either in a 
reactor or on an engine dynamometer yields the overall 
change in species concentration across the total length 
of the catalyst at discrete temperatures and flow rates for 
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a particular catalyst geometry and formulation providing 
solutions to Eq. 35, but only for a limited number of condi- 
tions. (It would be too time consuming and expensive to 
generate a catalyst map which would cover all possible 
operating conditions). We extrapolate the utility of the cat- 
alyst performance map by examining how the overall 
conversion is distributed locally along the catalyst length 
when a simplified reaction mechanism is considered to 
develop a scaling relationship between the catalyst map 
and the conditions being modeled. 

For example, suppose the reaction rate expression for 
each species i has the following simplified form: 



^(T w ,c: f ... > c^)=r r (c:)' 



(36) 



where r T is the reaction rate constant at 71, and n order is 



the reaction rate order. Substituting this expression into 
Eq. 35 and rearranging yields the following expression; 



(37) 

The numerical solution technique (presented later in this 
report) relies on the breakdown of the catalyst into a 
number of longitudinal computational elements. Conse- 
quently, it will be necessary to determine how the surface 
concentration changes on a cell by cell basis. This is 
obtained by integrating Eq. 37 over a the length of a sin- 
gle computational element to yield: 
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where Az y is the length of an individual computational 

cell at the j th location along the substrate, and C' w _ is the 
concentration of species t exiting cell j. 
Using Eq. 38 together with the definition of a local con- 
version efficiency, rfj t at the j th location along the sub- 
strate for species i, gives: 
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These expressions for local catalyst efficiency are related 
to the overall catalyst efficiency provided by the catalyst 

map, r( caJ , via the relationship shown in Eq. 40. 
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Note that in the case of w oftftfr =i, Eq. 39a yields a local 
conversion efficiency that is independent of inlet concen- 
tration and will be uniform along the entire active portion 
of the substrate when uniform longitudinal computational 
cell spacing is used (i.e., 

AZj -Az = L cat I N cell for all j ). For this special 

case, we define a uniform local efficiency as: 
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and rearrange Eq. 40 to give: 
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For a catalyst monolith with a square cell cross-section, 
Eq. 41 can be expressed as: 
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and the importance of the term in the exponential is that it 
gives the correct scaling relationship to extend our 
engine mapping data to catalysts of other geometries- 
other cell densities and wall thickness, differing space 
velocity and differing volume. This allows us to replace 
the detailed species conservation kinetic expression with 
a concise catalyst performance map. 
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One can also distribute the overall catalyst conversion 
efficiency for each species to each calculation node 

along the length of the catalyst by solving for j{uniform) 
and thereby estimate the energy release and species dis- 
tribution along the substrate. In general, ii } is a function 
of many variables. Here, we have chosen one functional 
dependence which corresponds to the fast mass transfer 
assumption where the overall conversion is controlled by 
a simple chemical reaction mechanism which in practice 
has produced results which compare favorably with 
experimental data obtained during typical emission cycle 
conditions. 

One can correlate other dependencies for special circum- 
stances such as high space velocity conditions, to simu- 
late an aging effect, or to mimic other unique features 
which occur during experimentation that are not well rep- 
resented by the above simplification. SIMTWC is written 
in a manner in which these can be easily integrated. 

With the above simplifications of fast mass transfer, local 
uniformity of conversion efficiency and quasi steady 
kinetics, the equation set reduces to 
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where we utilize our estimate of local efficiency to distrib- 
ute the exothermic chemical energy release discretely 
along the catalyst according to: 

fj Nred 
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where Q chem j is the exothermic heat released per 
"active" computational cell j (a computational cell is con- 
sidered "active" if it is above a specified light-off tempera- 
te): L catiaaive) ls the len 9 th of catalyst that is "active"; 
^ ceii(active) is the number of active computational cells 
assuming uniform spacing; N red is the number of reduc- 

tants in the calculation; and T w is the average surface 
temperature of the active cells. In SIMTWC, CO and H 2 

are considered as one species where C"- = X • C c g ° for 
gasoline [see Montreuil, Adamczyk, and Williams (1992) 
for a detailed discussion of why this is appropriate], and 
all HC are grouped as one species. In addition, to deter- 
mine a total oxidant level, NOx is lumped together with 
oxygen (with one mole of NOx equivalent to Vz mole 0 2 ). 
Thus, the energy liberated by reaction with 0 2 is consid- 
ered similar to that with NOx. Further, C 0 2, s tore refer s to 
two surface species, -- an oxidized site and a reduced 
site. They are related through a site conservation equa- 
tion specified as 



^0, Red _ ^Tot 



(50) 

This equation set is solved using the initial and boundary 
conditions shown earlier (Eqs. 22 - 27) that specify gas 
inlet temperature, flow rate, and species concentrations 
on a second-by-second basis during the specific drive 
cycle under consideration. 

The Oxygen Storage Model - As mentioned earlier, the 
oxygen storage mechanism in SIMTWC catalyst element 
is a transient phenomenon that removes oxygen from a 
lean exhaust gas stream and stores it on the surface of 
the catalyst. When the exhaust gas is rich relative to sto* 
ichiometry, the oxygen storage mechanism releases the 
oxygen stored on the catalyst surface to react with the 
reductants (CO and HCs) from the exhaust gas. The 
amount of oxygen that may be stored in the catalyst is 
called the total oxygen storage capacity and is a function 
of the amount of oxidizable material (e.g., ceria) in the 
washcoat of the catalyst. As the catalyst ages, the total 
oxygen storage capacity is reduced and the overall cata- 
lyst conversion efficiency generally declines with this 
reduction. 

In this model implementation, the rate at which the oxy- 
gen is stored and released is a function of the exhaust 
mass flow rate, the exhaust air-fuel ratio, the catalyst 
temperature, the oxygen diffusion between the gas and 
surface, the oxygen storage/release reaction rates on the 
surface of the catalyst, the oxygen storage capacity, and 
the amount of oxygen already stored on the catalyst sur- 
face. As seen above, the oxygen storage model is a first 
order model that accounts for these factors. Here, the 
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exhaust gas oxygen excess or deficiency relative to sto- 
ichiometry, 0 2icess , expressed as a mole fraction, is given 
by 



O 



(X-l) 



2xcess 



y/4 



+ 4.76 A 



(51) 



where y is the hydrogen-to-carbon ratio of the fuel, and 
X is the stoichiometrically normalized air-to-fuel ratio of 

the exhaust gas. A positive 0 2xccss (A > 1) indicates an 
excess of oxygen in the exhaust gas free stream while a 
negative 0 2xccss (A > 1) indicates a deficiency of oxygen. 

The normalized molar flow rate of excess/deficient oxy- 
gen into the catalyst, m , is given by 
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2 xcess ^exh 
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2, cap 
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where m exh is the exhaust molar flow rate and 0 l cap is 

the total oxygen storage capacity of the catalyst in moles/ 
Kg of catalyst. Using 



del 
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for each site species and the total 0 2 storage site conser- 
vation (Eq. 50) as the species conservation equation 
which describes 0 2 storage site action and defining 



Xo 2sbrt ~ 



C 



a 



2, cop 



(54) 



Then the overall rate of oxygen storage/release is given 
by 



d Xo Xm 
dt 



(/ = store &i = release) 



(55) 



where k i is a lumped parameter that includes the effects 
of oxygen diffusion from the gas to the catalyst substrate 
and the oxygen storage or adsorption/reaction (i=store) 
and release or desorption/reaction (i=release) rates; and 

Xo 2 ^ r is the fractional filling of the oxygen storage 
capacity (or oxygen bucket), 0 < Xo, - ' ■ When there 



is a deficiency of oxygen (A > 1), the'rate of change of 

Xo lyU „ is: 



dt 



X o->, Kt „ 



release 



(56) 



Thus, when the oxygen storage sites are nearly full 
( Xn ~ 1 ). * he oxygen release rate is near its maxi- 
mum. Conversely, when the oxygen storage sites are 
nearly empty (Xo, ~0), the oxygen release rate is 
near zero. Furthermore, when there is an excess of oxy- 
gen (A > 1) , the rate of change of (1- Xo, ut „ ) is propor- 
tional to (1 - ), 



\{ l -Xo^-Kore\^-Xo^ 



(57) 



or, 



dt 



Xo 2 „ ofr - K 



(58) 



Thus, when the oxygen storage is nearly full, the oxygen 
storage rate approaches zero. When the oxygen storage 
is nearly empty, the oxygen storage rate is near its maxi- 
mum. 

The effective relative air-fuel ratio, X cat , within the cata- 
lyst after the release or storage of oxygen remains to be 
calculated. The mole fraction of oxygen adsorbed or des- 

orbed by the catalyst in a time interval dt, 0 2 ads , is given 
by 



02.-*=- 



dt 



r 



WexH'dt 



If 



(59) 
When 



\0 2Md5 \>\0 Uccss \ then 0 2jub = 0 2m 
ads > 0 ' oxygen is stored by the catalyst and l ri 



becomes richer. Conversely, when 0, adi <0, oxygen is 
released from storage and the X cat becomes leaner. In 
both cases, X cat moves toward the stoichiometric value 
of 1 . The effective X cat is given by 
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(l-4.76(0 2xcess -0 2 ^)) 
When0 2 .^ = 0 2xccss , A Mf =l. 



(60) 



The Data Base Catalyst Map - The performance charac- 
teristics of the catalyst can be obtained either on an 
engine dynamometer or in a chemical reactor. These 
experiments provide overall conversion efficiency across 
an entire catalyst of prescribed geometry, PM formula- 
tion/loading and degree of aging. The steady state con- 
version efficiencies of CO, HC and NOx (defined as 

^cat '11 cat ^ca" . respectively) are correlated vs. a variety 
of inlet gas parameters that include temperature, air fuel 
ratio or lambda, flow rate, and emission gas concentra- 
tion level. Using the steady state data in conjunction with 

the model, the average surface temperature, T w , is 

calculated for each mapping data point and included in 
the map. In general, 



lcar = 



Icat (A>*T w9 m g ,C g , geometry, PM , aging , . . . ) 
where i equals CO, HC or NOx. 



(61) 



Ideally, one could experimentally determine these con- 
version efficiencies as function of all the above men- 
tioned parameters over a wide range of all possible 
conditions. However, to contain the total number of 
experiments and the corresponding facilities time to a 
practical value, the global conversion efficiencies of HC, 
CO and NOx are currently mapped for a single catalyst 
geometry and formulation as function of A/F ratio or 
lambda, for three inlet gas temperatures at three values 
of space velocity, and feedgas level concentration corre- 
sponding to a 4.6L 2 V Ford engine. 

Since the values of conversion efficiency are mapped 
over a finite range of conditions for a specific catalyst for- 
mulation, these efficiencies are interpolated when in the 
range of the experiments and extrapolated when outside 
the range at which they were obtained. The "core" part of 
a typical map is generated using inlet gas temperatures 

held constant at T map{min) = 380°C (720°F), T map{mi{i) = 

450°C (840°F) and T map{max) = 510°C (950°F), while con- 
tinuously varying A/F ratio from 1.0 A/F ratio lean of sto- 
ichiometry to 1.0 A/F ratio rich of stoichiometry at a 
space velocities of 85,000 hr"\ 170,000 hr"\ and 
225,000 hr* 1 . In addition to the core efficiencies obtained 
at constant inlet gas temperature, a "light-off" efficiency 
curve is generated which provides CO, HC and NOx effi- 
ciencies for a stoichiometric inlet gas which has its tem- 



perature slowly increased from ambient temperature to 
380°C (720°F). A complete map for a particular PM for- 
mulation also includes data obtained for several PM load- 
ings over the application range of these catalysts. 

In all mapping completed to date, the performance data 
were obtained after the catalyst sample was aged on the 
current Ford aging cycle for 125 hours of exposure. This 
level of exposure is equivalent to a certain level of vehicle 
driven miles which varies depending on the particular 
catalyst formulation. Use of the model to predict the per- 
formance of a catalyst at a different age requires genera- 
tion of another map. Nonetheless, a catalyst map for a 
single formulation and age can be used with appropriate 
scaling factors, discussed earlier, to account for differ- 
ences between the mapped catalyst and the modeled 
catalyst in terms of space velocity, and catalyst geometry 
(such as monolith cell density and wall thickness). 

Since our experiments yield global conversion efficiency 
for CO, HC and NOx and are obtained either at constant 
inlet gas temperature or at constant stoichiometry, we 
rewrite the global catalyst efficiency function in the follow- 
ing form which will lend itself to interpolation or extrapola- 
tion with a number of key variables. 



icai^f^M^SV.C'etc.) 



(62) 



where f scale is a scaling function (represented by Eq. 67 
in the next section) that is used to adjust the efficiency 

interpolated from the map (ijj^) for conditions different 

from those of the map. The map interpolation is per- 
formed in the following manner: 



n — 
'map 



J map * w / 



(63) 



where: 



fL P = core ma P efficiency 

fLp = Ugh* ~ °ff curve efficiency 

.map (rain) w w, map (max) 
~" -^vv. map (max) (f ^ ^v, map (max) 



(64) 



which means that we interpolate directly within the core 
experimental data base if T w - f . . < T < T m> (m lt 

^ w, map (mm) w w.mapivnjx) • 

and for surface temperatures outside the core experi- 
mental data base, an extrapolation is performed based 
on the shape of the ratio 

fmap(^>T w 

)! fLp^ T w.map(mm)) for surface tem- 



peratures lower 



map\ » w./nap(min) * 

than 



w,map( mint 



or 
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fLp^ T wj^t^) / fLpi^w^imBxi) for surface tem- 
peratures greater than T map{m ^ x) . This preserves the 
shape of conversion efficiency as a function of A/F ratio 
(or X) at all temperatures. // (A = 1,7^) is obtained 

from the catalyst lightoff efficiency data obtained at sto- 
ichiometry. Efficiency data are then scaled for other 
parameters such as space velocity, and inlet gas concen- 
tration. Finally, while the ideal catalyst map would include 
performance data for a variety of PM loadings and cata- 
lyst ages, these data are not always readily available. 
Consequently, to provide a first order estimate of the 
potential impacts of PM loading or age differences, the 
scaling function developed for SIMTWC also contains 
terms which can be used to provide some insight into 
how these differences can also affect catalyst perfor- 
mance. 

Scaling Function - The role of the scaling function is to 
extend the catalyst performance map to conditions occur- 
ring outside the performance map boundaries. This may 
include differences such as catalyst geometry, PM load- 
ing, age, inlet gas concentrations, and at the very least, 
flow rate (or space velocity). Since the catalyst map is 
obtained at a finite range of space velocities, a primary 
requirement of the scaling function is to predict perfor- 
mance at space velocities lower than 85,000 hr" 1 and 
greater than 225,000 hr* 1 . The following scaling for space 
velocity is applied and is based on the expressions spec- 
ified previously (see Eq. 42). 



(65) 



where rj 1 is the mapped catalyst conversion efficiency 

at the mapped space velocity SV ref , and Tj' cat is the esti- 
mated catalyst conversion efficiency for a normalized 
space velocity ( SV norm ): 



SV =3600 V @STP /(A r L ,SV f ) 

norm g \ face cat ref ' 



(66) 



We extend Eq. 65 to the following generalized form of the 
scaling function to cover a wide variety of possible differ- 
ences between the mapped catalyst and the catalyst 
being modeled. 



>L = i-(i-C) /sr 



(67) 

and SF' is the overall scaling factor for species i that has 
the following form: 



sr=sv norm 

(Age^.PM^.GSA^J 
(PPM^SVCR'J* 



(68) 



where: 



Age' Ave is the average reduction in catalyst performance 
due to aging (for each species i) for the entire catalyst 
based on a user supplied performance reduction profile 
for individual catalyst elements. 

PM' fac is an empirically based curve of the effect of PM 
loading and has the following format: 



(-cy PM / ) 
[ /™ load J 



(69) 



GSA accounts for differences in geometric surface 



area: 



(70) 



PPM' adj handles differences in inlet gas concentration 

that can affect catalyst performance, particularly in the 
case of species that follow a reaction mechanism more 
complex than first order. 



PPM 1 -(CI C V 



(71) 



and a final term, SVCR' is available to account for the 
varying effect space velocity has at different tempera- 
tures for each species i (included for future use, not cur- 
rently exercised) 

Available Oxidant Limitations - After the mapped effi- 
ciency is adjusted based on the scaling function, a final 
check is made to ensure that the reductant efficiencies 
have not exceeded the limit imposed by the amount of 
oxidant available. On the lean side, where there is an 
excess of oxygen, the reductant conversion efficiency is 
not affected by oxidant availability constraints. On the rich 
side, the deficiency of oxidants will place an overall limit 
on the possible conversion efficiencies of the reductants. 

Due to small errors in the mapping process, it is assumed 
in the model that a lack of oxidants to fully react with the 
reductants will primarily impact the CO conversion (HC is 
given priority to react with the available oxidant). The fol- 
lowing expressions are used to determine the maximum 
possible CO conversion imposed by oxidant limitations: 



9 



n co = ^-C; c (2 + y/2)77 

/ max 



HC 
cat 



C C0 + C H, 

Y = 

[(Cr + c;-')C + C(2 + y/2).r7-]. a 
a = l/Redox ratio /n rfd 
red,. 



Redox r „ ; „ = 



OX;. 



red in =C™ + C^ + (2 + y/2)-C l 

o*,,=2<=+c;°< 

r lrej=red bunud /red in 



HC 
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(72) 

(73) 
(74) 

(75) 

(76) 

(77) 
(78) 

(79) 



GENERAL CONSIDERATIONS -To analyze an emis- 
sions system using SIMTWC, component geometry and 
inlet gas temperature, A/F ratio, flow rate and species 
concentrations must be accurately specified. These 
feedgas values can be determined from other models 
which generate engine-out properties, or they can be 
obtained from vehicle experiments conducted over the 
appropriate drive cycle for analysis. Since all six inlet 
variables are necessary for calculation, the absence of 
an exhaust gas variable such as temperature will prevent 
accurate calculation. An average surface temperature 
and an average gas temperature are calculated in the 
catalyst element based on temperature distributions from 
the previous time step. Upon determining these average 
temperatures, an inlet value of A/F or Xj n | 6 t is used in cor- 
respondence with the 0 2 storage submodel to determine 
a new value for X within the catalyst after 0 2 is either 
added or removed from surface storage sites. After the 
new local value of X is calculated, the catalyst perfor- 
mance map is interrogated at the appropriate A/F ratio, 
surface temperature and space velocity and a global cat- 
alyst efficiency is determined. From the global efficien- 
cies, the local catalyst efficiency is determined based on 
the number of computational cells along the substrate 
and the species distribution is determined along the sub- 
strate. From this distribution of species, the local value of 

Qhemj is determined and the gas and wall energy bal- 
ances are calculated. 



Since the output from one pipe/catalyst combination can 
be used as input to a downstream pipe/catalyst combina- 
tion, a series of elements can be analyzed for any combi- 
nation of exhaust system components. This is typically 
the approach used during a complete analysis of an 
exhaust system configuration. 

EQUATION DISCRITIZATION - To solve the Eqs. 44-48, 
the length along the exhaust system is divided into num- 
ber of longitudinal computational elements and these 
equations are approximated by their corresponding finite 
difference relations. The representations of the partial 
derivatives are obtained from Taylor series expansions in 
space and time and these expressions are shown below. 



dt 



At 



dz Az 



dc' c^'-cr 1 



ci 



(80) 



(81) 



(82) 



For the conduction term in the solid, the second deriva- 
tive of T w with respect to z is written as: 



d 2 T (rj^C'+O 



dz 



-=B 

2 H t 



•0-/5,) 



(Az) 2 

(r -it +t" ) 

v w i-> w i *j,i ' 

(Az) 2 



(83) 

where p t is a weighting factor used to optimize calculation 
stability and accuracy. Normally, a value for 3 t of 0.5 is 
used during calculation [see Hauffmann, K.A. (1989)]. 
Substitution of the Taylor series expansion and weighting 
factors into Eq. 44 yields: 



„ +l _ +1) 



Az 



h > p K J/MC - C > + <» - P: xC - t:;[ )} 

J (84) 

where the change in gas temperature is equal to the 
weighted average heat flux from the gas to the wall. It 
should be noted that a spatial weighting factor, p z is also 
used to optimize the accuracy of the calculation. Similar 
substitutions into Eq. 45 yield: 
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yn+l nr-n 

P,A c p — 



[d- AK^ ■ j3,(r;;' -t; +i )]- 



Az 2 



(i-A)(C J , 1 -2r;+r; i ) 



(85) 



The boundary and initial conditions are discretized as fol- 
lows: 



Az 



/i+l 



y/i+l 

4z 



= 0 



(86) 

where N ce n is the last computational element along the 
substrate 



T;r=T g . inUI (t + At) (8?) 

" W (88) 

q(o,/)=c;, n( , ( (r) (8g) 

The above equations are applied at each computation 
node along the pipe or along the catalyst elements and 
the following system of equations is constructed to obtain 
a solution for T g and for T w . This system of equations is 
solved simultaneously to obtain T g and T w . The solution 
of this system is obtained using a Gaussian diagonaliza- 
tion and back substitution algorithm. 



J 8j-\ J J gj 

+D.T n+l + ET n * 1 + F.r* { = G 

J w j J t?j+\ J w j+l J 



(90) 



SOLUTION PROCESS FOR A SINGLE TIME STEP 

1. To start a calculation step, the gas, surface tempera- 
ture, and Q arrays are initialized using the values 
from the previous time step or as specified at the 
start of the simulation. 

2. The feedgas properties for a given time step of the 
simulation (temperature (K), A/F ratio, flowrate (kg/s), 
CO (kg/s), HC (kg/s), and NOx (kg/s) mass flow) are 
input, and the gas temperature for the first computa- 
tion node is set according to Eq. 87. 



3. An average surface temperature (T w ) (based on the 
temperature distribution from the previous time step) 
is calculated considering only elements above a 
specified light-off temperature using Eq. 49. If none 
of the elements is above the specified light-off tem- 
perature, T w is set equal to the maximum surface 
temperature within the catalyst. 

4. From the inlet value of A/F (or X), the exhaust gas 
oxygen excess or deficiency is calculated using Eq. 
51. 

5. Based on the result of Steps 3 and 4 and the degree 
of. oxygen storage present on the catalyst, the frac- 
tion of oxygen adsorbed or desorbed by the catalyst 
is calculated using Eq. 59. This value is combined 
with the result from Step 4 to determine a new local 

value of X cm within the catalyst using Eq. 60. 

6. A normalized space velocity (SV Mrm ) is calculated 
according to Eq. 66 based on the exhaust flow rate, 
catalyst volume, and the space velocity (SV r€ > f ) used 
during the generation of the catalyst performance 
map. 

7. The catalyst performance map is interrogated using 
the value of X cat calculated in Step 5, and the aver- 
age surface temperature from Step 3 to determine 

global catalyst map efficiencies (^ map ) for CO, HC, 
and NOx at SV ref using Eq. 63. 

8. The map efficiency, rj' map is checked and clipped (if 

necessary) to ensure that it is within a lower limit of 
zero and an upper limit of one. 

9. An overall scaling factor ( SF' ) is calculated for each 
species i using Eq. 68 that incorporates several indi- 
vidual effects into a single parameter for each spe- 
cies. These individual effects include: 

a. Normalized space velocity from Step 6 ( SV norm ). 

b. Catalyst age (Age' Ave ): average reduction in cat- 
alyst performance (for each species i) for the 
entire catalyst based on a user supplied perfor- 
mance reduction profile for individual catalyst 
elements. 

c. Precious metal loading effects on each species i 
(PM' fac ): calculated using Eq. 69 

d. GSA differences from the mapped catalyst 
(CSA ): calculated using Eq. 70 

e. Feedgas concentration levels for each species i 
( PPM' a(Ij )' calculated using Eq. 71 

f. The varying effect space velocity has at different 

temperatures for each species i (SVCR') 
(included for future use, not currently exercised) 
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10. The scaling function (Eq. 67) is used with the map 
efficiencies from Step 8 and the overall scaling factor 
from Step 9 to determine the overall global catalyst 
efficiencies possible without taking into consideration 
limitations imposed by the amount of oxidant avail- 
able. 

1 1 . The efficiencies determined in Step 10 are compared 
to the maximum efficiencies possible based on avail- 
able oxidant determined from Equation 72, and 
clipped if necessary to these maximum levels to pro- 
vide the overall global efficiency ( ri' cm ) for each spe- 
cies i. 

12. Cell-by-cell local cumulative efficiencies for each 
species i {rj' ieum ) up through cell j are calculated 

from the global efficiencies determined in Step 1 1 
using Eq. 91. 



<_-['-(«-*.r*-] 

where SV^ ective is determined from Eq. 92 



(91) 



SV* 

^ w effective 



J cum 

SF i /(2Age i J /Age i A J/j cl 



(92) 



and Age 1 -- reduction in cell j conversion perfor- 
mance for species i due to aging. 

13. The cell-by-cell local cumulative efficiencies for the 
reducing species, k, are used to calculate the exo- 
thermic heat, Qhemj * generated within each compu- 
tational element j using Eq. 93. 

4^,=|(nL-nL,)c;^(^) 

k 

forr]), cum <ll 

a*, m „=o, 

(93) 

where: AH k is the heat of reaction k. 

14. Knowing the energy released within each computa- 
tional element as determined from Eq. 93, the sys- 
tem of equations represented by Equation 90 is 
solved simultaneously to obtain the gas and surface 
temperature distribution along the catalyst elements. 
An iterative calculation is used to handle the effects 
of the temperature dependent gas and wall proper- 
ties (without updating the efficiency calculation or 
exothermic heat release) until a user specified con- 
vergence criteria is reached. 



INPUT DATA PHYSICS (FEEDGAS 
SPECIFICATION) 

A key input to SIMTWC is a precise characterization of 
the engine-out exhaust gas or "feedgas." More specifi- 
cally, SIMTWC requires feedgas mass flow rate, temper- 
ature, HC, CO, and NOx level, and A/F ratio (or oxygen 
content) as a function of time during a drive cycle of inter- 
est. Ideally, this information can be provided by either 
another simulation model, or actual test measurements. 
However, many engine feedgas simulation programs are 
valid only for fully warmed-up engine operation, and in 
many cases, it is the transient cold-start and warm-up 
portion of the emission cycle that dominates overall emis- 
sion performance. Therefore, in most cases, SIMTWC 
must rely on test measurements to provide the needed 
feedgas information over the drive cycle of interest. 

Several emission test facilities routinely measure most of 
the information required by SIMTWC. While the measure- 
ment techniques used to characterize the feedgas are 
sufficient for many types of emission assessment/calibra- 
tion work, there are several areas in which the measure- 
ments do not adequately reflect the true exhaust feedgas 
with the precision required by SIMTWC. Like a real cata- 
lyst, SIMTWC is extremely sensitive to minute changes in 
feedgas constituents, and in order to obtain meaningful 
results from SIMTWC, it is very important that the input to 
the model reflect the actual feedgas quality as closely as 
possible. 

By understanding what causes the inaccuracies in the 
feedgas measurements, it is possible to make some 
adjustments to the reported feedgas values to more 
accurately reflect the true exhaust make-up. Five known 
causes of measurement inaccuracies are discussed 
below together with techniques to compensate for them. 

TEMPERATURE MEASUREMENTS - Thermocouples 
are the typical temperature sensors used to measure 

exhaust gas temperature ( T g ) in dynamometer and vehi- 
cle experiments. They are comprised of the junction of 
two wires of dissimilar materials which contains thermal 
mass. The temperature reported by the thermocouple is 

actually that of the of the thermocouple junction ( T TC ) 
and not necessarily that of the surrounding environment 
due primarily to the thermal inertia of the junction mass. 
The response lag of the thermocouple to changes in the 
exhaust gas temperature is a function of the thermocou- 
ple size and exhaust flow rate. By using a simple model 
of thermocouple energy conservation as suggested in 
Figure 2 f it is possible to account for the effects of ther- 
mocouple mass and to obtain a better estimate of the 
actual exhaust gas temperature. This is done by taking 
into account the heat transfer between the exhaust flow 
and the thermocouple junction and writing the corre- 
sponding energy balance between the thermocouple 
junction and the gas flow. This yields: 
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dT TC _ 
dt it^cCPtc 



i T rc-T g ) 



(94a) 



1 = hjAn; 

dT, 



T g = T TC + r 



TC 



dt 



(94b) 

where x, in general, is a function of the exhaust flow rate 
and thermocouple size. 

Thermocouple Energy Balance 

Flow 




Figure 2. Schematic of thermocouple energy exchange 
process. 

Experience has shown that for the range of exhaust flow 
rates encountered during a standard EPA Federal Test 
Procedure (FTP), x is nearly independent of flow rate, 
and use of a constant value for a given thermocouple 
size is sufficient for calculation. Values for x range from 4 
seconds for a 1/16" thermocouple to 8 seconds for a 1/8" 
thermocouple. For simulations of cycles that involve 
higher exhaust flow rates (such as the US06 in the Sup- 
plemental FTP), it is necessary to take into account the 
effect of flow rate on x. 

EXHAUST FLOW MEASUREMENTS - Exhaust flow 
measurements at many emission test facilities are made 
using the C0 2 tracer method. This method relies on a 
measurement of engine out C0 2 concentration obtained 
from an analyzer which is located remotely from the test 
vehicle (-10 m). Due to the inherent time delay associ- 
ated with the transport of the exhaust sample to the emis- 
sion analyzer, the C0 2 tracer technique generally 
underestimates the engine flow dramatically during the 
first 5 to. 15 seconds of an emission test. This underesti- 
mate of flow is extremely important, since it also leads to 
an underestimate of the calculated energy transport from 
the engine to the emissions system, thus lengthening the 
calculated time of catalyst lightoff in the emissions sys- 
tem. It is therefore necessary to correct the exhaust flow 



determined with the C0 2 tracer technique during the 
early part of a test using other exhaust measurement 
techniques, such as the Smooth Approach Orifice (or 
SAO) or other engine parameters such as rpm, load, and 
manifold vacuum - a speed-density calculation of flow. 

Additionally, in many facilities, flow estimates from the 
C0 2 tracer technique are reported as the dry exhaust 
flow rate (total exhaust flow minus the water content). 
Because typical gasoline exhaust contains about 12-13% 
water, use of the dry exhaust flow rate would result in a 
significant underestimation of the amount of sensible 
heat delivered to the catalyst during the light-off portion 
of the emission test. Consequently, the dry exhaust flow 
is increased to account for the water content before the 
data is fed into SIMTWC. 

STOICH I0METR IC AIR/FUEL RATIO 

DETERMINATION - Transient Air/Fuel ratio (A/F) is gen- 
erally calculated from data obtained using conventional 
emissions analyzers or measured with a UEGO (univer- 
sal exhaust gas oxygen) sensor during an emission test. 
It is not unusual for measurements using these devices to 
be shifted slightly rich or lean (± .05) of the exhaust's 
actual A/F as a result of inaccuracies in instrument cali- 
bration, fuel hydrogen/carbon fluctuations, sampling 
delays, and only partial mixing of the exhaust gas from all 
engine cylinders. Inaccuracies of this order may appear 
small compared to typical exhaust A/F of 1 4.6:1 for gaso- 
line. However, when attempting to determine the relative 
oxygen content of the exhaust which involves the differ- 
ence between the measured A/F and the stoichiometric 
A/F, measurement shifts of this order can have a tremen- 
dous impact on the accuracy of calculation. Conse- 
quently, it is necessary to carefully analyze each data set 
to determine whether there is an offset between the mea- 
sured A/F ratio and the actual A/F ratio. 

In our process, an initial estimate of the offset in A/F ratio 
is obtained by averaging the measured A/F ratio values 
during a section of the test when the engine is known to 
be fully warmed-up and remains in closed-loop operation. 
At this condition, the A/F ratio should be near to the sto- 
ichiometric A/F ratio as determined from the H/C and O/C 
ratios of the fuel. If there is a difference between the mea- 
sured average A/F ratio and the stoichiometric A/F ratio 
calculated from fuel properties, the difference is elimi- 
nated for the entire emissions test. In addition, if a base- 
line emission test with tailpipe measurements is 
available, it is possible to make refinements to this initial 
estimate by comparing the ratio of tailpipe CO/NOx mass 
predicted by SIMTWC to that actually measured. 

TIME ALIGNMENT OF SIGNAL TRACES - During an 
emission test, the exhaust sample passes through sev- 
eral different sensors/analyzers to determine the mass 
flow, temperature, and A/F, HC, CO, and NOx levels. 
Generally, an attempt is made at the emission facility to 
time align the signals from the various detectors; how- 
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ever, in most cases, a more refined alignment is required 
to get the most meaningful results from emissions mod- 
els. Currently, we take advantage of the strong correla- 
tion between exhaust mass flow, temperature, and 
feedgas NOx mass flow rate that generally exist to align 
these three signals. We then assume that the HC and CO 
signals are aligned with the NOx, and finally use the 
strong correlation between CO concentration and A/F to 
align the A/F with all signals. 

"FILTERED" EMISSION SIGNALS -In many emission 
facilities, the CO, HC, and NOx concentrations are mea- 
sured by analyzers which are located at least 10m away 
from the test vehicle. Because of the long sample lines, a 
degree of mixing or filtering of emissions constituents 
occurs by the time the exhaust sample reaches the emis- 
sion analyzers. In experiments using a closely-located 
HC analyzer, known as the "Fast- FID" [Summers T. and 
Codings, N (1995)], this instrument has shown repeatedly 
the dramatic difference between the actual exhaust HC 
concentration and that measured by the standard remote 
analyzer. Again, to get the best possible results from 
SIMTWC or other emissions models, it is imperative to 
know the emission concentrations entering the catalyst 
as precisely as possible. 

Use of the Fast-FID is helpful, but it is not standard prac- 
tice at most facilities. In addition, other "Fast" instruments 
are needed. Fast-NOx analyzers are another example of 
instrumentation which shows dramatic differences with 
conventional analyzers and are being evaluated. In the 
interim, we are investigating methods to "unfilter" the 
emission signals using techniques similar to that used for 
the thermocouple signal. However, at this time, a stan- 
dard technique has not yet been adopted. 

RESULTS AND DISCUSSION 

MODEL COMPARISON WITH VEHICLE 
EXPERIMENTS -To assess modeling accuracy, 
SIMTWC was used to calculate catalyst performance and 
the corresponding accumulated emissions over several 
vehicle test cycles for a variety of different test vehicles 
and catalyst configurations as listed in Table I. In addition 
to those listed in Table I, other variable changes included 
catalyst position, catalyst cell density, catalyst PM load- 
ing, exhaust system geometry, catalyst aging, and fuel 
type (California reformulated and Indolene). All calcula- 
tions were done with fixed calibration of the model to 
determine the amount of difference between experiment 
and model calculation over a wide range of conditions 
and emissions levels. 



Table I. Vehicle and Drive Cycle Combinations 



Vehicle 


Engine 


Drive Cycle 


Grand Marquis 


4.6L V8 


FTP 


Town Car 


4.6 LV8 


FTP, MVEuro, & US06 


Contour 


2.0L 14 


FTP 


Heavy-Duty 
F-Series Truck 


6.8LV10 


Heavy-Duty 
Gas Engine Test 


Lincoln LS 


4.0L V8 


FTP 



Figure 3 shows the values predicted by SIMTWC when 
compared to those measured during experiment. The 
values on the abscissa represent the accumulation of 
transient experimental "modal" data obtained on a sec- 
ond-by-second basis throughout a typical vehicle-on 
chassis dynamometer experiment. The ordinate repre- 
sents the accumulated tail pipe emissions as calculated 
by SIMTWC (open square markers). The solid straight 
line would represent perfect agreement between experi- 
mental results and theoretical predictions. The two outer 
dotted lines indicate a difference of 30% between the 
modal measurements and the SIMTWC prediction. 
Herein, most predictions by SIMTWC lie within 30% of 
the measured values for all combinations of test condi- 
tions and geometries that generate emission levels that 
span several orders of magnitude in value. 

To help put this into perspective, Figure 3 also shows the 
corresponding "bag" emission values measured during 
these experiments (solid triangle markers). Note that the 
values from the bag measurements vary from the values 
of the measured modal results by approximately the 
same degree as the model prediction, thus suggesting 
that the error generated from the model is similar to the 
error which is generated during experiment. Moreover, 
remember that these bag measurements were obtained 
during the same experiment as the accumulated modal 
measurements, using the same experimental facility and 
vehicle emissions system, and still the difference in mea- 
surements is similar to the difference between model and 
experiment. Also, remember that the catalyst representa- 
tion in the model reflects the performance of an average 
aged catalyst and does not exactly represent the catalyst 
on any one particular vehicle due to manufacturing vari- 
ability, differences in catalyst aging (ADP/fleet), degree of 
contamination, and so on. Performance differences due 
to differences in aging are also discussed later. 
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Figure 3. 
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Figure 4 shows two vehicle geometries, which are deriva- 
tives of the Lincoln TownCar production geometry, used 
during calculations to compare second-by-second CO, 
HC and NOx emissions at catalyst midbed and tailpipe 
positions. 

Using configuration A and the appropriate feedgas input, 
the results in Figure 5 show the accumulated CO, HC, 
and NOx emission level during bags 1 and 2 of FTP cycle 
as a function of time at both midbed and tailpipe loca- 
tions. These accumulated sec-by-sec emissions show 
good correlation between the predicted and measured 
values at the midbed catalyst location. Here, they are 
within 30% of one another, with the difference in mea- 
sured and calculated HC within 10%. However, at the 
tailpipe position, the measured NOx emission is 50% 
lower than the calculated NOx emission. This is in part 
due to the fact that the lead catalyst was aged in the 
same manner as the catalyst used to generate the perfor- 
mance map in the model, whereas the second and third 
catalysts on the vehicle were aged behind the first cata- 
lyst as a system and experienced significantly different 
aging conditions than the first catalyst (e.g., lower inlet 
gas temperatures, little exotherm on the catalyst and 
reduced chemical poisons). This combination of factors 
suggests less deterioration of the rear catalyst bricks and 
is consistent with the higher NOx emissions produced 
during calculation. Since the current model calibration 
yields higher NOx emissions for rear catalyst position (a 
lower NOx efficiency), it acts as a conservative estimate 
of the tailpipe emission for the system. In addition, this 
suggests that an additional calibration should be under- 
taken for catalyst in the rear brick position and should be 
dependent on the temperature history of these catalyst 
which, in turn, is dependent on their location in the emis- 
sions system. 
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Figure 4. Exhaust system geometry 
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Accumulated emissions over the FTP drive 
cycle at midbed and tailpipe positions for 
TownCar Config. A. 



In Figure 6, the instantaneous CO, HC and NOx conver- 
sion efficiencies ({[]jn - DoutV Bin x 100) are plotted versus 
time at the midbed catalyst position for the first 200 sec- 
onds of the FTP drive cycle for TownCar configuration A. 
In general, they follow similar trends with experiment. 
However, for the first 10 to 20 seconds of the cycle, the 
measured efficiencies rose to unreasonably high values 
considering the catalysts were cold when the experiment 
started, and they (the efficiencies) dropped to negative 
values (not plotted) within a short period of time after 
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vehicle start. This behavior suggests the absorption and 
release of unreacted emissions species by the catalyst 
during its initial warm up. Since these experiments show 
little occurrence of reaction, the absorption/desorption 
simply acts as a time delay of emissions species without 
change in their total value and this effect is currently not 
considered in the SIMTWC model. A future SIMTWC ele- 
ment will be developed to account for this phenomenon. 
Many of the differences incurred between measured and 
calculated efficiency after this initial warm up period are 
due to the lack of precise alignment of AF ratio with emis- 
sion measurements and are also due to the effect of fil- 
tered measurements of emission signals. Moreover, 
when the feedgas concentration of a given emissions 
species is small, the error in calculating efficiency is 
greater. This is partly due to the accurate measurement 
of a small emissions concentration and its placement in 
the denominator during the efficiency calculation. As 
noted in the previous plot this has little effect on the accu- 
mulated mass emission from the vehicle, since the 
feedgas emission levels are extremely low during these 
periods. 
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Figure 6. Instantaneous CO, HC and NOx efficiency at 
midbed and tailpipe positions for TownCar 
Config A. 

As in Figure 5, Figures 7 and 8 show a comparison 
between calculated and measured accumulated CO, HC 
and NOx midbed and tailpipe emissions when drive cycle 
is changed (Figure7; TownCar configuration A; US06 
drive cycle) and when emissions system geometry (Fig- 
ure 8; TownCar configuration B: FTP drive cycle) is 
changed. As seen, the absolute emissions levels have 



changed in correspondence with driving cycle and with 
geometry. 

In Figure 7, the calculated NOx emission is in good 
agreement with measurement on the US06 drive cycle 
which does not contain a vehicle cold start. The calcu- 
lated HC emissions is within 30% ({1 - .68/.95}x100) of 
measurement. Since the regulation for the US06 drive 
cycle is a composite of the HC and NOx emissions and 
since the NOx emissions, as seen in Figure 7, is gener- 
ally an order of magnitude greater than the HC emission, 
the calculated US06 HC+NOx composite emission is in 
excellent agreement with measurement. 
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Figure 7. Accumulated emissions over the US06 drive 
cycle at midbed and tailpipe positions for 
TownCar Config. A. 

Figure 8 shows a comparison between the calculated 
and measured emissions generated by TownCar configu- 
ration B. Here, geometry has changed from the baseline 
configuration A. In this configuration, the catalyst is closer 
to the engine, the volume of the catalyst is greater (46 vs. 
42 in 3 ), its cell density is higher (470 vs. 400 cpsi), the 
catalyst shape has changed (round vs; oval), the feedgas 
signatures are different, and it has been aged at a differ- 
ent time than the baseline configuration. Yet, SIMTWC 
has accounted for these changes. As in Figure 5, Figure 
8 shows smaller differences across the lead catalyst and 
larger differences across the rear catalyst. Again, this 
catalyst experienced significantly different aging condi- 
tions than those used to age the lead catalyst and those 
used to age the catalyst to generate the performance 
map for the model. As seen in Figure 8, the largest differ- 
ence is in the NOx tailpipe emission; it is within 35% ( {1 - 
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.08/.12}x100) of measurement. Again, the current cali- 
bration in SIMTWC yields a conservative (slightly lower 
NOx efficiency) estimate of the tailpipe emission. 
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8. Accumulated emissions over the FTP drive 
cycle at midbed and tailpipe positions for 
TownCar Config. B. 



MODEL SENSITIVITY - SIMTWC contains all the sensi- 
tivities of catalyst and emission system behavior. As an 
example, if one biases A/F ratio by an extremely small 
amount, the emissions system efficiency will change 
appropriately. If geometry, or materials of the catalyst 
change, the appropriate material properties are incorpo- 
rated into the model. If engine spark timing or air/fuel 
ratio is varied as a method of engine control, the model 
reflects this engine control, since the feed gas time signa- 
ture used as input to SIMTWC has incorporated the 
appropriate information to account for this control. In this 
section, we give a few examples of performance changes 
with catalyst location change and A/F ratio bias change. 

Using TownCar configuration B as a baseline and its cor- 
responding feed gas data set over bags 1&2 of the FTP 
drive cycle, the changes in tailpipe emissions levels and 
emission system temperature are shown in Figures 9-11 
as catalyst location is changed and as A/F ratio is biased 
relative to stoichiometry. 

Figure 9 shows the effect of changing the position of the 
lead catalyst from its base line location (1 .5" downstream 
of the manifold flange, see Figure 4) to 0.3m ahead of the 
manifold flange and to 0.6m downstream of it. Note that, 
when the position of the first catalyst was moved down- 
stream and reached the position of the second catalyst, 
the entire catalyst system was relocated downstream 



appropriately. As seen in Figure 9, there is a 6x (0.06/ 
0.01) increase in the tailpipe HC emission, a 2x (0.7./ 
0.35) increase in the CO emission, and a 50% ({0.16/ 
0.11 -1}x100) increase in the NOx emission when the 
lead catalyst location is changed over this range. The 
large increase in the HC and CO emissions is due to 
retarded catalyst lightoff which is due to less sensible 
heat arriving at the catalyst. The smaller increase in NOx 
is due to lower temperature at the catalyst during steady 
state operation. 
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Figure 9. Emissions with change in catalysts 
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Figure 10. Feedgas temperature into 1st catalyst used in 
Figure 9. Higher temperatures correspond to 
the shortest distance from the engine. 
Distances correspond to those used in Figure 9 

Figure 10 shows how the feedgas temperature to the 
lead catalyst changes when the location of catalyst is 
changed as above. As seen, when the lead catalyst is at 
its closest position, the inlet gas temperature is approxi- 
mately 15% higher than when it is at its most down- 
stream location. This is extremely important, since it 
affects catalyst lightoff, hence emissions, and also affects 
catalyst deterioration due to high temperature exposure. 
Therefore, a tradeoff between lightoff performance and 
catalyst deterioration ultimately determines the optimum 
location for the lead catalyst, since the emissions system 
must perform well at 120K vehicle miles. 
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Since a small bias in A/F ratio strongly affects emissions 
performance, we show the effect of AF ratio bias on 
tailpipe emissions for this catalyst configuration. Here, we 
artificially shift A/F ratio either lean or rich of the baseline 
A/F ratio indicated as zero bias on Figure 1 1 . As seen in 
Figure 11, the CO emission increases 3x (1.2/0.4) when 
A/F is biased rich by 0.02 A/F and the NOx increases by 
80% ({0.18/0.1-0.1}x100) when the A/F ratio is biased 
lean by a corresponding amount. Remembering that 
these increases are specifically related to the exact A/F 
ratio-time signature studied, the amount of increase or 
decrease can change dramatically from these levels. An 
important point to remember is that the specification of A/ 
F ratio by experiment or by other upstream models of 
engine performance or control strategy must be within 
±0.1% (0.02/14.6) of the actual A/F ratio or the tailpipe 
emissions can be in error by 100%. This accuracy sug- 
gests the importance of A/F ratio specification or mea- 
surement as all calibration specialists know. 
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Figure 1 1 . Effects of Air-to-Fuel ratio bias on tailpipe 
emissions over the FTP drive cycle for 
TownCar Configuration A. SUM indicates the 
sum of all three emissions constituents after 
normalization with the ULEV emissions 
targets (HC=0.04; CO=1.7; NOx=0.2). 

In addition to the variation of each individual emissions 
species, the sum of all three species after normalization 
with the ULEV emissions targets (HC=0.04g/mi; 
CO=1 .7g/mi; NOx=0.2g/mi) is shown. Thus 

„„, CO HC NOx 

SUM = + + 

1.7 0.04 0.2 

and is non-dimensional. The rescaling of emissions spe- 
cies is done to provide equal weighting before summa- 
tion. This value shows a minimum very near the zero bias 
point, suggesting that the overall A/F ratio measurement 
is correct and the vehicle calibration is appropriate for 
minimum emissions over the FTP cycle. 



CONCLUDING REMARKS 

A simple tailpipe emission prediction model (SIMTWC) 
has been developed by combining fundamental conser- 
vation expressions for mass and energy with a first order 
oxygen storage model and an empirical t data base of 
steady-state catalyst performance obtained during 
engine dynamometer assessments. SIMTWC offers sev- 
eral advantages over more complicated techniques 
including computational expedience and robustness with 
accuracy usually better than 30% for all emission species 
and over various drive cycles. Using this approach, an 
FTP emission cycle for a six pipe plus six brick emission 
system can be modeled on a 233 MHz Pentium™ II PC 
in less than three minutes. This provides an engineer the 
ability to rapidly evaluate many geometries and operating 
strategies. Updates to SIMTWC as catalyst formulations 
change are also accomplished in a relative easy and 
timely manner making SIMTWC is an extremely useful 
tool which offers the ability to sort through major after- 
treatment design and strategy alternatives quickly while 
assisting in the interpretation of experimental data. 

ACKNOWLEDGMENTS 

The authors thank Mr. J.D. Pakko and Mr. A. Kolasa for 
obtaining the catalyst performance mapping data on the 
engine dynamometer. They thank Ms. A. Bork for con- 
ducting reactor experiments and Mr. D. Osborn for oper- 
ating the dynamometer cell. The authors also thank Mr 
R.J. Pawlowicz and Mr. S.C. Williams for supplying vehi- 
cle data for comparisons with model results. 

REFERENCES 

1. Gandhi, H.S., Delosh, R.G., Piken, A.G., and Shelef, 
M. t "Laboratory evaluation of three-way catalysts," 
S.A.E. Transactions, Vol 85, Sec. 2, p. 201, The Soci- 
ety of Automotive Engineers, (1976). 

2. Koltsakis, G.C., Konstantinidis, PA., Stamatelos, 
A.M., "Development and application of mathematical 
models for 3-way catalytic converters," Applied Catal- 
ysis B: Environmental, Vol 12, p 161-191, (1997). 

3. Hauffmann, K.A., Computational Fluid Dynamics for 
Engineers, Engineering Education System(TM), Aus- 
tin, TX, 1989. 

4. Li, P., Adamczyk, A.A. and Pakko, J.D., "Thermal 
Management of Automotive Emission Systems: 
Reducing the Environmental Impact", Thermal Engi- 
neering for Global Environmental Protection, ISBN 1- 
56700-055-X, Begell House, Inc, Publishers, New 
York, Wallington (UK), p 55-77 (1996), Ed. S. Sen- 
gupta. 

5. Montreuil, C.N., Williams, S.C. and Adamczyk, A. A., 
"Modeling current generation catalytic converters: 
Laboratory experiments and kinetic parameter opti- 
mization -- Steady state kinetics," S.A.E. Paper No. 
920096, The Society of Automotive Engineers, 
(1992). 



18 



6. Oh, S.H., Cavendish, J.C., and Hegedus. L.L t "Math- 
ematical modeling of catalyst converter lightoff: sin- 
gle pellet studies/ AiChE Journal, Vol 26, No. 6, p 
935, (1980). 

7. Otto, N.C. and LeGray, W.J., tt Mathematical models 
for catalytic converter performance," S.A.E. paper 
No. 800841, The Society of Automotive Engineers, 
(1980). 

8. Summers T., Collings, N., 'Modeling the Transit Time 
of a Fast Response Flame Ionization Detector During 
In-Cylinder Sampling," S.A.E. Paper No. 950160, 
The Society of Automotive Engineers, (1995). 

NOMENCLATURE 



internal heat transfer coefficient augmenta- 
tion factor (dimensionless) 
external heat transfer coefficient augmenta- 
tion factor (dimensionless) 

facial area of the catalyst substrate (m 2 ) 

Hgt CC + Tt Hgt 2 1 4.0 

(racetrack geometry) 

thermocouple surface area (m 2 ) 

cross-sectional area of pipe wall (m 2 ) 



face 

K - 

w iat 

CY 

^ 1 PA4 



C2\ 



PM 



cy 



c 



ppm 



c[ 



ci 



g. inlet 



c 



store " 
->Tot 

^O-, store ~ 



C 



C 

»V inlet 



c 

W. outlet 

cc 



= solid cross-sectional area of catalyst (m 2 ) 

average reduction in catalyst performance 
due to aging for species i 
precious metal loading function constant for 
species i 

precious metal loading function constant for 
species i 

precious metal loading function constant for 
species i 

initial catalyst concentration of species i at 
time t=0 

constant for species i used in the inlet con- 
centration adjustment factor 
gas concentration of species I 
(g-mole/m 3 @ STP) 
inlet gas concentration of species I 
(g-mole/m 3 @ STP) 

0 2 storage concentration of species i 
(g-mole/m 3 @ STP) 

Total 0 2 storage concentration (g-mole/m 3 
@ STP) 

surface concentration of species i at location 
j (g-mole/m 3 @ STP) 

surface concentration of species i at the inlet 
of the catalyst (g-mole/m 3 @ STP) 
surface concentration of species i at the out- 
let of the catalyst (g-mole/m 3 @ STP) 

W - Hgt (m) (racetrack geometry) 



CD 
CPrc 

c Pw 

C P „ icer) 

'cat 
•pipe 

f 

J map 
f 

J Tmp 
f scale 

8 

Ga 

Ga map 

Gr 

CSA,.„„ 
Hgt 



cell density (cells/m 2 ) 

heat capacity of the exhaust gas 
(J/Kg/K) 

962.097 + 0.1507 -T g ( K) 
thermocouple heat capacity (J/Kg/K) 

wall heat capacity (J/Kg/K) 
heat capacity of ceramic (J/Kg/K) 

3.435- 10 7 



1071.0+0.156- T W (K)-- 



TJK) 2 



kL 



= heat capacity of iron (J/Kg/K) 

= 460.52 

= heat capacity of stainless steel (J/Kg/K) 

= 502.39 

= heat capacity of steel (J/Kg/K) 

= 502.39 

= catalyst cell hydraulic diameter (m) 

- ^K"^) 

= inside pipe diameter (m) 

= outside pipe diameter (m) 

= core map efficiency of species i 

= light-off curve efficiency of species i 

scaling function to extrapolate conversion 
efficiencies outside performance conditions 

= gravitational acceleration (m/s 2 ) 

= 9.807 

= geometric surface area (m 2 /m 3 ) 

= 4.o(i.o-w; 4cd)4cd 

geometric surface area of mapped catalyst 
" (m 2 /m 3 ) 

= Grashof Number 

ratio of modeled catalyst geometric surface 
area to mapped catalyst surface area 
height of substrate face (m) 
(racetrack geometry) 

= internal heat transfer coefficient (J/m 2 /s/K) 

heat transfer coefficient between flow and 
thermocouple (J/m 2 /s/K) 

= external heat transfer coefficient (J/m 2 /s/K) 

= release rate constant for oxygen storage 

= storage rate constant for oxygen storage 

= mass transfer coefficient of species i (m/s) 
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length of catalyst (m) 

length of catalyst that has a temperature 
above a specified light-off value (m) 

length of pipe (m) 

normalized molar flow rate of excess/defi- 
cient oxygen into the catalyst 

exhaust gas molar flow rate (g-mole/s) 
exhaust gas mass flow rate (Kg/s) 
inlet exhaust gas mass flow rate (Kg/s) 

thermocouple mass (Kg) 

exhaust gas molecular weight (Kg/g-mole) 
0.028966 

reaction rate order 

number of computational cells per exhaust 
system component 

number of computational cells having a tem- 
perature above a specified light-off value 

number of reactions 
number of reductant species 

number of species 
Nusselt Number (inside) 

Laminar Nusselt Number (inside) 
Turbulent Nusselt Number (inside) 
Nusselt Number (outside) 

measure of total oxidant concentration within 
the catalyst (g-mole/m 3 @ STP) 
oxygen storage capacity (g-moles 0 2 /Kg cat- 
alyst) 

mole fraction of oxygen adsorbed or des- 
orbed by the catalyst 
exhaust gas oxygen excess or deficiency 
(mole fraction) 

Gas pressure (Pascal) 

total catalyst cell inside perimeter (m) 

4 O,,., CD A /°c< 

catalyst outside perimeter (m) 
it Hgt + 2 CC (racetrack geometry) 
pipe inside perimeter (m) 

' ft* 

pipe outside perimeter (m) 



PMi 



foe 
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red in 
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Re 

Redox„,„ 



sv 



° v effective 



SV 

° V norm 
SKef 

SVCR { 
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Hz 



empirical function which describes the effect 
of PM loading on species i conversion 

= PM loading (Kg/m 3 ) 

= concentration effect factor for species i 

= Prandtl number of exhaust gas at T g 

rate of heat generated by chemical reaction 
per unit length of catalyst (J/m-s) 
rate of heat generated in cell j by chemical 
reaction per unit length of catalyst (J/m-s) 

measure of total reductant concentration 
consumed (g-mole/m 3 @ STP) 
measure of total reductant concentration 
within the catalyst (g-mole/m 3 @ STP) 

- reaction rate constant at T w 

= gas constant 
= 8.314 J/g-mole/K 

= Reynolds number 

rate of reaction of the k th reaction per unit 
surface area (g-mole/m 2 /s) 
overall rate of reaction for reactions involving 
species i (g-mole/m 2 /s) 
= redox ratio 

= red in lox in 

= external surface-to-volume ratio (m 2 /m 3 ) 

= P oJ^ 

overall scaling factor for species i 
(dimensionless) 

= space velocity (hr -1 ) 



= 3600- V, 



@STP 



'(A fc 



face ^cat 



K D, 



effective space velocity taking into account 
aging and the overall scaling factor (hr* 1 ) 
normalized space velocity (dimensionless) 

sv/sv ref 

reference space velocity (hr* 1 ) used during 
the mapping of the catalyst, 
space velocity vs. temperature adjustment 
factor for species i 
film temperature (K) 
(T g + T w )/2 

exhaust gas temperature at position z and 
timet(K) 

^Vin/«(0 = ,n,et 9 as temperature at time t (K) 
T mQp = Inlet mapping gas temperature (K) 

T rc = thermocouple junction temperature (K) 

T , _ surface temperature at position z and time t 

" (K) 
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Initial wall temperature at position z (K) 
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Initial temperature distribution at time=0 
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average catalyst surface temperature of cells 
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calculated average catalyst surface tempera- 
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ture under mapping conditions (K) 
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ambient temDerature (K\ 


0 / r v 

r w(Ft) 


- 


t 
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local gas velocity along the catalyst (m/s) 
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volumetric gas flow rate (m 3 /s) 
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width of substrate face (m) - 
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(racetrack geometry) 
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cell wall thickness (m) 
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fuel hydrogen to carbon ratio 
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axial position along length of the pipe or 
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catalyst (m) 
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normalized position along length of the 
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catalyst (dimensionless) 
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Greek 
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a 




available oxidant calculation variable 




= 






expansion coefficient (1/K) 
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P 




reciprocal of absolute gas temperature 


T 








1/T amb 






P, 




temporal weighting factor 











spatial weighting factor 
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void volume fraction 
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(l.O-W, yfCDj 
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local efficiency of species i in cell j 


A T 1 • 

\H HC 
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nL 




overall catalyst efficiency for species i 










local efficiency for species i in cell j based on 






T/( inform) 


rfi umfo^ = constant for all cells 






n' 

'map 




performance map efficiency for species i 






rt' 




maximum conversion efficiency for species i 






1 max 




based on available oxidant 






1r e J 




total reductant conversion efficiency 






Y 




available oxidant calculation variable 







- exhaust gas density (Kg/m 3 ) 



PMw p 

g g_ 

RT„ 



= 352.55/ T g (K) 



wall density (Kg/m 3 ) 

density of ceramic substrate (Kg/m 3 ) 
1985 

density of Iron (Kg/m 3 ) 
7593 

density of metallic substrate (Kg/m 3 ) 
8121.1 

stoichiometrically normalized air-to-fuel ratio 
stoichiometrically normalized air-to-fuel ratio 
within the catalyst after adjusting for oxygen 
storage 

thermal conductivity of gas (J/m/s/K) 
8.459x1 0 _3 +5.7x10" 5 T g (K) 

wall thermal conductivity (J/m/s/K) 

thermal conductivity of stainless steel 
(J/m/s/K) 

11.416+.0130T W (K) 
thermal conductivity of metallic substrate 
(J/m/s/K) 
31.0- ,0095T W (K) 

thermal conductivity of Iron (J/m/s/K) 
56.653- .01 732T W (K) 
thermal conductivity of ceramic substrate 
(J/m/s/K) 

0.0506(32.734 - .01 T W (K)] 
dynamic viscosity of exhaust gas (N-s/m 2 ) 
[.00268T g (K)+1.384]*10- 5 
thermocouple time constant (seconds) 
fractional filling of the oxygen storage 
capacity 

heat of reaction per mole for the k th reaction 
(J/g-mole) 

heat of reaction for CO (J/g-mole) 

2.83 x 10 5 

heat of reaction for H 2 (J/g-mole) 

2.84 x 10 5 

heat of reaction for HC (J/g-mole) 
6.37 x 10 5 

= uniform computational cell length (m) 
= j lh computational cell length (m) 
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